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ABSTRACT 

In this paper, an economical low-noise plasma simulation model origi- 
nated by Denavit is applied to a series of problems associated with electro- 
static wave propagation in a one -dimensional, collisionless, Maxwellian 
plasma, in the absence of magnetic field. In Part I, the model is described 
and tested, first in the absence of an applied signal, and then with a small 
amplitude perturbation. These tests serve to establish the low-noise fea- 
tures of the model, and to verify the theoretical linear dispersion rela- 
tion at wave energy levels as low as lo"*"' of the plasma thermal energy; 
better quantitative results are obtained, for comparable computing time, 
than can be obtained by conventional particle simulation models, or direct 
solution of the Vlasov equation. The method is then used to study propaga- 
tion of an essentially monochromatic plane wave. Results on amplitude 
oscillation and nonlinear frequency shift are compared with available 
theories. The additional phenomena of sideband instability and satellite 
growth, stimulated by large amplitude wave propagation and the resulting 
particle trapping, are described in Part II. 
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1 . INTRODUCTION 


This paper is concerned with computer simulation of electron plasma 
waves in a Maxwellian plasma, in the absence of magnetic field. The aims 
of our simulations are, first, to develop an economical low-noise simula- 
tion technique, and second, to apply it to the study of linear and non- 
linear wave phenomena in a one-dimensional plasma with periodic boundary 
conditions. Throughout the work, emphasis is placed on simulations for 
wave and plasma parameters comparable with those assumed in available 
theories, and accessible in laboratory experiments. 

There are two distinctly different approaches to the simulation of 
plasma dynamics: first is the use of a particle simulation model in 
which individual charged particles are followed, and second is direct 
numerical solution of the Vlasov equation describing the charged parti- 
cle velocity distribution function. The particle simulation model has 
the disadvantage that the fluctuation level is usually several orders of 
magnitude higher than in an actual plasma. This stems from the fact that 
it is not feasible to follow on the computer the dynamics of as many parti- 
cles as there are in a plasma. The fluctuations not only give rise to non- 
physical effects, but also make it difficult to study linear and weakly 
nonlinear phenomena. This is particularly unfortunate since most of the 
nonlinear theories to date are based on an expansion method which is valid 
only in weakly nonlinear cases . Consequently, they cannot be clearly vali- 
dated by particle simulation, nor vice versa . Direct solution of the 
Vlasov equation is subject to numerical instability associated with the 
free- streaming term in the Vlasov equation. This tends to limit application 
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of the method to short-time simulations. However, since the Vlasov equation 
does not contain discrete particle encounters and thermal fluctuations of 
macroscopic quantities such as the electric field and charged particle 
density, the behavior of very small amplitude waves can be studied in a 
quantitatively accurate manner, 

Denavit has proposed a hybrid approach which can drastically reduce 
the fluctuations inherent in the particle simulation model, and greatly 
ease the computational difficulties of the Vlasov equation approach. His 
model is a particle simulation model in the sense that the motions of a 
large number of particles are followed in time. The particles do not keep 
their identities, however, and there are similarities to the Vlasov approach 
in that values of the velocity distribution function are defined on a grid 
in phase-space. Despite its attractive features, and obvious potential for 
application to a wide range of linear and nonlinear problems, little use of 
the approach seems to have been made so far. The purpose of this paper is 
to apply it to a logical series of such problems, involving electron plasma 
waves in a Maxwellian plasma, and so provide results for quantitative com- 
parison with available theories and experimental results. 

In Section 2, a hybrid simulation model suitable for studying one- 
dimensional electron plasma waves is described. In Section its equili- 
brium characteristics are tested, i.e., noise growth with zero applied 
signal amplitude is examined. In Section k, linear wave propagation is 
studied, for a very small amplitude monochromatic signal. Progressive 
increase in signal amplitude causes amplitude oscillation, nonlinear 
Landau damping, and nonlinear frequency shift effects to appear. These 
phenomena are studied in Section 5 . The results are discussed, brief ly 
in Section 6. 
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2. HYBRID SIMULATION MODEL 


The model to be used in this work employs a cloud-in-cell (CIC) scheme. 
This was originally devised to reduce fluctuations in particle simulation 
models, and involves considering the particles as clouds of distributed 
charge; the center of the cloud is taken to be the particle coordinate, 
and a spatial grid is used for computing field quantities. 2 In addition 
to this grid, we introduce a grid in velocity-space, and represent the 
particles by points in (x-v) phase-space as shown in Fig. 1. The phase- 
space is consequently covered with a rectangular grid of dimensions 
, Av . The velocity grid extends from v to v , where v and 
V 2 are chosen such that the numbers of particles with velocities in the 
intervals v < v 1 or v > v 0 are negligible. 

In creating a plasma with a Maxwellian velocity distribution, our 
model employs the following method. The particles are equally divided 
into a number of velocity groups. All the particles in one group are 
assumed to have the same velocity, v , and mass and charge are assigned 
to them in proportion to exp (-v^/2v 2 ) , where v fc is the electron 
thermal velocity. This is shown schematically in Fig. 1. Since the 
charge-to-mass ratio is the same for all of the particles, the accelera- 
tion is also the same. One of the advantages of this method of generating 
a Maxwellian distribution by weighted particles is that improved resolu- 
tion is provided in the tail of the velocity distribution, compared with 
a Maxwellian distribution with identical particles. 

The system is set up at time t = 0 using a quiet start technique, 
and proceeds as in a CIC model. After a certain number of time-steps. 
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the distribution function, is reconstructed at the grid points by periodic 
smoothing, and is interpreted as representing a distribution of new dis- 
crete particles. The motions of these particles are followed until the 
next reconstruction. The quiet start technique, and the periodic smoothing, 
are essential parts of the hybrid approach proposed by Denavit. They are 
used to achieve a very low fluctuation level, and allow the model to be 
applied to a wide range of linear and nonlinear problems. 


2.1 Quiet Start 

The quiet start technique, proposed by Byers, is a method of elimi- 
nating macroscopic fluctuations in a particle code at early stages of 
evolution.^ This is done by placing the particles only on the equilibrium 
trajectories in phase-space at time t = 0 , as shown in Fig. 1, and in 
principle provides a noiseless system. In practice, round-off errors due 
to the finite number of digits representing the numbers in the computer 
introduce some fluctuations. However, this level is many orders of magni- 
tude lower than that in the particle codes . 

Although the quiet start technique works well at early times, ft 
ceases to be effective after a time, 2n7(k Av) , where k is the ma xi - 
mum wavenumber possible in the system. 1 This breakdown occurs because the 
velocity distribution of the particles is being replaced by a set of dis- 
crete beams. Such a system is also subject to streaming instability, even 
if the envelope of the beam density is Maxwellian.^ In the limit of small 
beam spacing, the temporal growth rate, ox , is given by 1 


CD. 

X 




( 1 ) 
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Therefore, a simulation can be carried out with very low fluctuations up to 
-1 

times of order ox . Since ox-*-0 as Av->0 , in principle it is possible 
to perform a low-noise simulation for as long as is desired by choosing 
a sufficiently small Av . In practice, however, Av cannot always be 
made as small as is desirable, because the smaller Av is to be, the 
more particles are necessary. Periodic smoothing may be used to combat 
this instability, and to make a long-time simulation possible with a 
relatively small number of beams.* - 


2 .2 Periodic Smoothing 

Periodic smoothing constitutes a periodic averaging of the distri- 
bution function in phase-space. It can be expressed by 



v ' )dx'dv' 


( 2 ) 


where w and w are weighting functions for coordinate- and velocity- 

X V 

r*j 

space, f is the averaged distribution function, and the integration is 
over the whole of phase-space. 

In particle models with a phase-space grid such as that shown in 
Fig. 1, the integral in Eq. (2) reduces to a sum over the collection of 
particles, and we want to find the averaged distribution function, f , 
at the phase-space grid points. If f(x',v%t) is taken to be the mass 
of a finite-size particle, the center of which is located at (x',v') at 
time t , then Eq. (2) implies that the value of f at the (i-j) grid 
point, (x.,Vj) , is obtained by distributing the mass of each particle 

among the neighboring grid points according to the weighting prescribed 
by w x and w v . This is a reconstruction of the distribution function 
from a given distribution of particles. 
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Since the smoothing process expressed by Eq. (2) does not conserve 
all. of the moments of the distribution function, it causes diffusion of 
the distribution function. It is to be expected that, with judiciously 
chosen weighting functions, this diffusion may suppress the streaming 
instability with a minimum of other undesirable effects. Such weighting 
functions are described, and their diffusion rates are estimated, in 
Reference 1. In our work, we shall use the quadratic weighting func- 
tion which conserves particles, momentum, and energy. The frequency of 
smoothing necessary to quench the streaming instability depends on the 
choice of weighting function and the beam spacing. A rough estimate 
would be that at least one smoothing operation is necessary within the 
growth time of a perturbation with the maximum wavenumber of the system. 

2 .J Comparison with Denavit's Model 

In the work of Denavit, 1 the Lewis variational method'’ was used to 
construct a model. In the Denavit model, finite-size particles are 
chosen to have a triangular spatial distribution, instead of the uniform 
charge distribution of the CIC model. The numerical scheme based on that 
model turned out to be more complicated, therefore more time-consuming, 
than our model. In addition, although his scheme is energy-conserving, 
it does not conserve momentum, whereas the CIC scheme is formulated in 

g 

such a way as to conserve momentum. The non-conservation of momentum 
indicates the existence of self-force, i.e., a particle is effected by 
the force due to the field that is created by itself, which is non- 
physical. The energy-conserving feature of the Denavit scheme may not 
be very useful in practice, since energy conservation is exact only in 
the limit of a vanishingly small time-step. 
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It may be remarked, further, that in testing his model Denavit chose 
to study two-stream instability. He compared his results with those from 
a particle code 7 9 and the Vlasov approach. 9 These simulations were car- 
ried out at relatively high electrostatic energy levels, i.e., lo"^ - lo " 2 
times the total energy, which is to be compared with an order of 10~^ in 
our tests to be described in Section 4 . 
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3 ■ EQUILIBRIUM BEHAVIOR 


It is logical to begin our series of numerical experiments on a 
Maxwellian plasma with no applied perturbation at all. Here, and 
throughout the paper, only electron motions are followed; positive ions 
are considered to constitute a homogeneous, immobile, neutralizing back- 
ground. Periodic boundary conditions are applied in space: a particle 

leaving one end of the system is immediately reintroduced at the other 
with the same velocity. The most important parameters are the number 
of time-steps, N g , after which the smoothing operation is repeated, 
and the beam spacing, Av 

In Eig . 2, the total field energy is plotted against time for 
Av = v t /7 . It will be seen that by increasing the frequency of 
smoothing, i.e., decreasing N g , the streaming instability is sup- 
pressed; for < 32 , the field energy stays roughly constant 

throughout the simulation run. To study this further, the initial 
energy spectrum, and the time-averaged energy spectrum are shown in 
Fig. 3. The energy spectrum for N g = 32 shows that mode energy 
tends to increase as time increases. For N g = 1 6 and 8 , the mode 

energy seems to stay at roughly the same level, i.e., the streaming 
instability is stabilized. The total energy of the system, i.e., the 
particle kinetic energy plus the field energy, was found to be con- 
served to within 0.1* up to cn p t ~ 270 , where m p is the electron 
plasma frequency. 
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Time-averaged energy spectrum of the system shown in 
Fig. 2. 
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In a simulation for Av = v t /lh , it was observed that the growth 

rate of the field energy is greatly reduced compared with the case 

Av = v t /7 > and larger values of N are adequate to suppress the 
l s 

streaming instability 

The important fact established here is that, for stable cases, the 

O 

total field energy is fluctuating at a level 10 times lower than the 

thermal energy of the plasma during the whole run. Since the level of 

fluctuations in the particle model for a system of length L = J2 X^ , 

where X^ is the electronic Debye length, with N = k0$)6 electrons, is 

-5 -2 

of the order of 10 — 10 times the thermal energy, we have achieved 

50 - 60 dB reduction in fluctuations. The fluctuation amplitude in the 
particle model can be reduced by increasing N , but it should be remem- 
bered that the level varies only as N -1 / 2 . 
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4. LINEAR WAVE PROPAGATION 


The purpose of this simulation was to verify predictions of Landau 
damping for electron plasma waves. Waves were excited at t = 0 by 
applying the perturbations 

Ax i = A\ d cos kx i , Av^ = Av fc sin kx^^ , ( 5 ) 

where x^^ , Ax^ , and Av,^ are the position, displacement, and velocity 
perturbation of particle i , A is the amplitude, and k[= STm/L] is 
the wavenumber. In this simulation, only Mode 2 (n = 2, i.e., there are 
two wavelengths in the system) was excited. The results are shown in 
Fig. 4. For N < 16 , there is excellent agreement with the theoreti- 

cal prediction by Langdon,^ shown by solid lines, which takes into account 
finite-size particle effects and spatial grid effects. 

Although the fluctuation amplitude due to round-off errors in Fig. b 
increases with time, note that the ratio of electrostatic energy to ther- 
mal energy, [(eE/m^) /v fc ] , at t = 0 is 4.25 X 10 in this simu- 
lation. Particle simulation with such good quality, at such a low electro- 

static energy level, has not been feasible with previous models. For 
example, in order to reduce the fluctuation level to 10 times the thermal 

energy in a particle simulation with the same system length, it would re- 

3 4 

quire 10 - 10 times more particles than are used in this simulation. 

Since the computing cost increases roughly in proportion to the number of 
particles, it would be prohibitively expensive. In contrast, the computing 
cost in this simulation was found to be less than twice that with the corre- 
sponding CIC model. Suppose the smoothing operation is performed every 
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FIG. U. Simulation of Landau damping. Solid lines are the prediction of 
the Langdon theory Dashed line is the prediction of point 
particle theory. 
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16 time-steps. One smoothing operation in our computer code takes about 
7 sec on an IBM 360/6 7 for 8 I 92 particles. It is equivalent to an increase 
of about 0.44 sec per time-step. Since it takes our computer code about 
0.75 sec per time-step for the CIC model with 8 I 92 particles and a system 
with 128 cells , the total computing time per time-step is about 1.2 sec. 

In addition to our check on temporal Landau damping of a signal , we 
have verified the predicted linear dispersion characteristics of electron 
plasma waves . The results are shown in Fig. 5, and agree well with theory. 
The initial perturbations were applied for Modes 7-20 with random phases. 

The electrostatic energy of the individual modes excited was about 

-6 

4 X 10 times the thermal energy at t - 0 

The simulation results presented in this section serve to demonstrate 
that the hybrid approach provides quantitatively accurate results on the 
collective behavior of plasma in the linear regime, where comparison with 
theory can readily be made. There is no reason why it should not be an 
equally effective tool in the nonlinear regime, for which analytical 
results are not so readily available. In the remainder of the paper, we 
shall employ it in the study of a number of nonlinear wave phenomena which 
appear when the signal amplitude is increased. Before doing so, however, 
we may comment on a well-known phenomenon associated with simulation by 
the Vlasov equation approach. 

4.1 Recurrence Phenomenon 

In a simulation such as that carried out in this section, a perturba- 
tion with wavenumber k first damps to a low level at the Landau damping 
rate, and then reappears suddenly at time t = 27r/kAv with higher 
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amplitude than its initial one. After its reappearance, the perturbation 
decays with a damping rate slightly different from the Landau damping rate. 
This recurrence results from approximating the continuous distribution 
function by a finite number of delta-function beams: the perturbation with 
wavenumber k on each beam comes back into phase after a time t r = 27r/kAv , 
and with larger amplitude than its initial value due to the streaming insta- 
bility. In the limit of infinitely many beams, phase-mixing prevents recur- 
rence of the initial state. The recurrence phenomenon was observed by 
Denavit.^ Similar phenomena have been observed in the numerical solution 

of the Vlasov equation by the Fourier-Hermite method, 1 ' 1 ' the finite dif- 
11 12 

ference method, and the Lewis variational method, and have been ascribed 
to finite resolution in describing velocity-space oscillations . According 
to our own numerical experiments, periodic smoothing alone does not suffice 
to prevent the recurrence. 

In anticipation of succeeding sections on nonlinear phenomena, we can 
say that the recurrence phenomenon has not been found to pose any problems 
in our simulations. In Section 5, on large amplitude wave propagation, no 
irregularity is demonstrated at the recurrence time, t r . Similarly, 
where growth of small amplitude test waves is concerned, no irregularity 
is observed. 
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5 . NONLINEAR WAVE PROPAGATION 


Since plasma is a highly nonlinear medium, the linearized analysis 
gives only a limited description of its behavior. The question arises 
of what will happen to Landau damping and wave dispersion when the wave 
amplitude is increased. Theoretical studies of this question were first 

V 

"t ’2? 1 1 

made by O'Neil, ^ and Al'tshul and Karpman. These authors found that 
the amplitude changes in time in an oscillatory manner, after initial 
Landau damping. The amplitude oscillation is due to periodic exchange 
of energy between the wave and electrons trapped in the potential wells 


of the wave. The exchange occurs on a time scale of l/o^ , where 

1 /2 

(ekE Q /m e ) ' ] is the bounce frequency of an electron oscillating 


at the bottom of a potential well, k is the wavenumber, and 
the wave electric field amplitude. 


is 


In solving the Vlasov equation, O'Neil, and Al'tshul and Karpman, 
assumed that the amplitude variation is so small as to satisfy 

« 1 > where is the Landau damping rate. Under the same 

assumption, Bailey and Denavit have taken into account the effects of 


the slowly-varying amplitude, and obtained essentially the same ampli- 
tude oscillation, except that the time at which the amplitude begins 
to grow again after the initial damping is delayed. ^ Gary has treated 
the case 7^/^ > ^ analytically, and shown that the wave starts to 
decay at a rate smaller than the Landau damping rate at a time when the 
linear theory is expected to break down. 1 ^ Sugihara and Kamimura have 
considered a wide range of TjV 0 ^ values. Unlike the theories mentioned 
so far, their treatment is self-consistent: it includes the interaction 
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between the electric field and the averaged electron velocity distribution. 

Recent work by Oei and Swanson is also self-consistent , and gives similar 

18 

results to those of Sugihara and Kamimura for 0 < < 1 

All of the theoretical studies discussed above assume that the elec- 
tric field is so small that the distribution function in the resonant 
region can be expressed by a Taylor expansion about the wave phase velocity 

up to the first order term in velocity. This condition may be written as 

2 

» ( V p/ V t ) ; where v^ is the phase velocity of the wave (o*/k) 

This is such a stringent condition that it is not easy to meet in laboratory 

experiments, especially when also satisfying the condition « 1 

A few experimental data on Landau damping of large amplitude waves 

19 

have been reported. Malmberg and Wharton observed spatial amplitude 

1^5 

oscillation in qualitative agreement with the O’Neil theory y modified to 

20 

fit the spatial case by Lee and Schmidt . Oei and Swanson compared their 

theoretical results with the experiments of Malmberg and Wharton, and found 

agreement on the amplitude oscillation lengths but not on the detailed 

18 

behavior of the amplitude. One of the reasons may be that their experi- 

2 

mental parameters do not meet the condition au/o^ » ( V p/ V t) - Speci- 

2 

fically, they have y^/o^ ~ 0.1 and co/o^ ~ ( V p/ V t ) for th e results which 
exhibit amplitude oscillation. Franklin ejt a_l. have made detailed measure- 
ments of the spatial dependence of amplitude for electron plasma waves with 

2i 

different initial amplitudes, i.e., for different values of . How- 

ever, for large initial amplitude, they failed to obtain results in agree- 
ment with the theory. This was ascribed to the appearance of sideband 

22 

growth due to trapped particle instability.""^ Their experimental parameters 
for the measured results, corresponding to < 0.^5, yield < l;(v h 
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This suggests that comparison of the available theories with the experiments 
is inappropriate. 

In view of the foregoing difficulties, computer simulation suggests 
itself as a means of bridging the gap between the theoretical assumptions 
and readily attainable experimental parameters, It allows conditions to be 
studied for which analytical approaches are not tractable. Such simulations 

23 2b 

have been carried out by Knorr, " and Armstrong, using direct solution of 

25 

the Vlasov equation, and by Dawson and Shanny, using the particle simula- 
tion model. Knorr observed a decrease in the damping rate for large ampli- 
tude waves at times such that ~ 1 . Armstrong considered the same 

problem and found in addition to Knorr T s results that large amplitude waves 
grow again after damping initially. He also found that the initial damping 
of a large amplitude wave is stronger than is predicted by the Landau theory. 
A similar observation of the enhanced initial damping was made by Dawson and 
Shanny . 

One of our purposes has been to use the hybrid model described in 
Section 2 to investigate the nonlinear behavior of longitudinal monochromatic 
plasma waves more comprehensively than has been possible previously. Another 
has been to investigate the nonlinear frequency shift of electron plasma 
waves. In a plasma of infinite extent, or of finite length with periodic 
boundary conditions, the frequency of a wave of large amplitude deviates 
from that of a small amplitude wave due to nonlinear effects. In an experi- 
mental plasma, in which a wave is excited at a fixed frequency, the shift 
should occur in wavelength instead of frequency. ' 

The frequency shift has been studied analytically by Manheimer and 

2 6 27 28 2Q 

Flynn, Morales and O'Neil, Dewar, and Lee and Pocobelli, y and found 
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to be proportional to E ' ' . So far, there has been no report of labora- 

tory observations of nonlinear wavelength shift for comparison with these 
theories. In this section, we shall test the theoretical predictions of 
nonlinear frequency shift against computer simulations carried out by use 
of the hybrid model . 


5*1 Computations 

We have performed a series of computer simulations to demonstrate the 


nonlinear behavior of monochromatic electron plasma waves in a collision- 
less plasma. The electrostatic energy of the waves in these simulations 
was of the order of 10 4 times the thermal energy. This is about two orders 
of magnitude smaller than in the simulations of Dawson and Shanny. 2 ^ Some 

23 2k 

of the simulations by Knorr,^ and Armstrong, are in our range of energy. 
Their computations have not, however, been carried out for long enough times 
to demonstrate amplitude oscillation. 


Amplitude Oscillation ; Figure 6 demonstrates clearly this phenomenon. 
Mode 3 was excited initially according to Eq. ( 3 ), and the evolution of the 
amplitude was followed in time with periodic boundary conditions applied in 
space. The amplitude oscillates and approaches a constant value due to 


phase-mixing of the trapped particles, and formation of a Bernstein-Green- 
Kruskal (BGK) mode."^ 


The initial amplitude of the wave was eE„/m v 0) ~ 5.4 X 10~ 2 

O' e t p - > 

corresponding to a bounce frequency of o^/o> ~ 0.09 • The measured 

initial damping rate is 7 L /^ p S O.OII 9 . These combine to give 
7 L Afe 2 0.13 . The measured frequency is ~ 1.15 . The corre- 


sponding wave phase velocity is v /v fc ~ 3 . 9 I , so that oi/o^ ~ 13 and 

(v p /v fc ) ~ 15 are of the same order. 
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FIG. 6 . Amplitude oscillation of a large amplitude monochromatic wave. [Av = v /l4, N = 16 

t/ * g 

N = 8192, L = 6k Ax, H = Ax = Vj. = -4.25 v t , v 2 = h.Qs At = 0 . 25 /co ] . 
(Nonsymmetric velocity-space is used to provide resonant particles at high velocities.) 




^-t&^re 7 shows the temporal behavior of the distribution function, in 

the vicinity of the wave phase velocity , averaged over space. The changes 

in the distribution are relatively small; for example, at oj t = 48 the 

P * 

ratio of the peak value to the maximum value of the (nearly) Maxwellian 
distribution is of the order of 10 A bump is formed in the distribution 
for ui p t ~ 48 , and reappears for u^t ~ 144 and 2kO . Comparison with 
Fig. 6 indicates that these times correspond approximately to minima in 
the amplitude. The height of the bump becomes progressively smaller on 
its reappearances, because of phase-mixing of the trapped particles. ^ 

A similar bump was observed by Armstrong, and considered to cause growth 
of waves with phase velocities lying in that region of the bump that has 

24 

positive slope. The bump on the tail of the distribution function has 
spatial structure. This may be contrasted with the initially spatially 
homogeneous distribution whose evolution is considered in the quasilinear 
theory of a warm beam-plasma system.^ 1 

In Fig. 8 , we present the results of a series of simulations for various 
values of the initial electric field, , expressed in terms of the con- 
venient parameter 7 L /^ , where we recall that = (ekE^inJ 1 / 2 . Only 

one mode was excited at t = 0 for each simulation run, and a different 
mode and amplitude were used in each run. The amplitude was normalized to 
unity at t = 0 in the plots, it will be seen from Fig. 8 that amplitude 
oscillation occurs for small values of 7 ^/ 0 ^ , and that the oscillation 

becomes less pronounced, with Landau damping extended for a longer period, 
as 7 l /o^ increases. The fluctuations in the curves for large values of 
are ^ ue to round-off errors made in representing numbers by a 
finite number of digits in the computer. 
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FIG. 7. Temporal behavior of the spatially-averaged distribution function in the simulation 
shown in Fig. 6. The phase velocity of the wave is marked by an arrow. 
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FIG. 8. Temporal evolution of amplitude for various initial amplitudes, i.e., 



Frequency Shift : Figure 9 shows the variation of the nonlinear fre- 

quency shift as a function of the electric field amplitude. Mode 5 was 

excited initially according to Eq. ( 3 ), with amplitude (eE_/m v m ) 

O' e t p ' 

varying from small values (9 X 10 which exhibit Landau damping, to 

large values (3 X 10 ) such as were studied in the simulations of Dawson 

25 

and Shanny. For each simulation with different amplitude, the frequency 

of Mode 3 was measured by computing the total amount of phase change in 

the Fourier transform of the electric field between 00 t = 6 and 60 . 

P 

The frequency shift plotted in Fig. 9 was then obtained by subtracting 
the linear frequency, = 1 .21+7 , obtained from the Langdon theory , 0 

from the measured frequency. Except for very small amplitudes, the non- 
linear frequency shift is proportional to , and given by 
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To check the dependence of this result on the beam spacing, the simulations 
were repeated with the beam spacing halved, and the same number of smooth- 
ing operations. The differences in frequency shift were not more than 3 $. 

A significant fact to note here is the high degree of accuracy with 
which it was possible to determine the frequency, and frequency shift. The 
model based on the hybrid approach is, therefore, much more efficient than 
a particle code in terms of computing cost for this measurement. 


5 -2 Comparison with Theory 

First, we may check the observed initial damping against the linear 

theory of Langdon, ^ which includes finite size particle and spatial grid 

effects. The theoretical values of the Landau damping rate, y and 

L ^ 
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Nonlinear frequency shift of an electron plasma wave 
L = 50 X_, H = Ax = (50/6U)X_, v, = - 3.79 v , v n = 5 




frequency, , obtained by retaining only the p = 0 term in the summa- 

tion in the expression for plasma permittivity given in Reference 6 are 
7 ^/^p = 0.0118 and o^/a^ = l.lV? for Mode J plotted in Fig. 6 . We see 
very good agreement with the measurements described in Section 5.1. The 
theoretical predictions for each mode presented in Fig. 8 have also been 
found to agree with the measured initial damping and frequency with errors 
of less than 1 $. 

Amplitude Oscillation ; Next we consider the nonlinear theories due 

13 13 17 

to O’Neil, Bailey and Denavit, and Sugihara and Kamimura . O'Neil 
obtained a theoretical expression for the time-varying damping rate by 
using an energy conservation relation between the resonant particles and 
the wave. Substituting the initial damping rate ( 7 ^/ 0 ^ = 0.0119) 
obtained from our computer simulation (Fig. 6), and the bounce frequency 
(o^/oip = 0 . 09 ) calculated from the initial amplitude in the same simulation, 
we obtain the amplitude variation plotted in Fig. 6 . After damping initially, 
the wave starts to grow somewhat earlier than it does in the simulation. This 
can be ascribed to the variation in wave amplitude, which was not taken into 
account in calculating particle trajectories. 

Bailey and Denavit incorporated the effects of slowly-varying wave ampli- 
tude to lowest order in d/a 2 , where a(t) = [ekE(t )/inJ = «( 0 ) , 
and a. ~ da/dt . We have solved numerically a set of differential equa- 
tions which they obtained, for the same values of 7 ^ and 0 ^ used above, 
and with the results plotted in Fig. 6. There is very good agreement between 
the theory and the simulation. We note, however, that there is a slight dif- 
ference in amplitude, and that the phase— mixing is somewhat slower in the 
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simulation results than the theory predicts. These differences are prob- 
ably due, first, to the fact that the condition, cd/o^ » (v 
not satisfied in the simulation, and second, that the theory of Bailey and 
Denavit is not self-consistent. 


P /v t r 


Sugihara and Kamimura derived from the Vlasov equation a set of integro- 
differential equations which describe the behavior of the amplitude of a mono- 
chromatic wave. Numerical solutions of these integro-dif f erential equations 
demonstrated amplitude oscillation for y^/u^ « 1 > and Landau damping 

for 7 l A^ » 1 • Some of their results are reproduced in Fig. 8. First, 

we note that their calculation shows that, for =0.1 , the amplitude 

approaches a constant value after nearly two periods of oscillation, although 
the distribution function still retains nonuniform features. In our simula- 
tion, however, the amplitude oscillation lasts more than two periods, and 
does not seem to die out so quickly. This fact seems to be in at least 

qualitative agreement with a nonlinear spatial Landau damping experiment by 

19 

Malmberg and Wharton 7 in which there was no clear sign of phase-mixing. A 
similar feature of this persistent amplitude oscillation was also observed 
in the behavior of an externally excited large amplitude wave in a simula- 
tion of sideband instability by Denavit and Kruer.^ Second, we recall that 
Sugihara and Kamimura found that there is a critical value of y = 0.77 } 
which separates waves into those with oscillatory behavior (7 L A^ < 0.77) , 
and those which are continuously damped > 0.77) . Figure 8 indi- 
cates that there is no such critical value below y^/o^ = O.93 . Third, 

we note that there is a tendency in our simulation results for the amplitude 
to decrease to a lower level, for a given value of , than is predicted 

by the theory of Sugihara and Kamimura; the first maximum is also lower than 
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the theory predicts. Although the simulation results given in Fig. 8 are 
similar to the theoretical results obtained by Sugihara and Karaimura, it 

p 

is important to note that in our simulations > (v /v t ) , whereas 

P 

they implicitly assumed that cd/o^ » (v /v fc ) 
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Frequency Shift ; Manheimer and Flynn examined the self-consistency 

J -z 

of the O’Neil solution for the time-asymptotic state: J they studied 
whether the potential created by the 0 'Neil solution satisfies the Poisson 
equation. They found that it is approximately self-consistent if a frequency 
shift given by 



is included, where g is a numerical factor equal to 2*^, f is the 
initial distribution function, and is the linear plasma permittivity. 

In deriving Eq. (5), Manheimer and Flynn only considered the trapped par- 
ticles with simple harmonic motions, i.e., those near, the potential wells 
of the wave, and the untrapped particles with straight line orbits. Morales 

and O'Neil solved an initial value problem to find the time -dependent shift 

27 

in the complex frequency of the wave. They took into account the exact 
trajectories for both the trapped and untrapped particles, and obtained a 
frequency shift which varies in an oscillatory manner and approaches a 
constant value in the time-asymptotic limit. Their time-asymptotic frequency 
shift is expressed in the same form as Eq . (5) except that g «s I.63 


Lee and Pocobelli predicted frequency shifts for waves with v /v > 4 

P t ~ 

up to about 50$ larger than those predicted by Morales and O'Neil. These 
were obtained by including effects of electrons not in the vicinity of the 
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phase velocity of the wave ." 7 In contrast to these theories treating the 
case in which the wave is switched on suddenly at t = 0 , Dewar considered 

the case of an adiabatically excited wave, i.e., the wave was turned on 

28 

gradually. He obtained a time-asymptotic frequency shift similar to that 
expressed by Eq. (5), but with 0=1.09 • 

Substituting <\h v = 1.2V7 for Mode 3 , obtained from the Langdon 
theory and the Maxwellian distribution for f in Eq . ( 5 ), we have 

-0.19^ (Morales and O'Neil), 

(6) 

-O.ljo^ (Dewar) , 

which are plotted in Fig, 9* We see that the slopes of the lines from the 
simulation, and from the theory of Morales and O'Neil, are very similar. 

This is to be expected because our simulation of an initial value problem 
resembles the Morales and O'Neil problem, rather than the Dewar problem. 

It should be remembered, however, that the theoretical result is the time- 
asymptotic value, whereas the measured frequency shift is an average over 
the period u^t = 6 to 60 . It should also be recalled that the value of 

corresponds to the initial amplitude of the wave. Since the theoreti- 
cal result due to Morales and 0 'Neil was obtained under the condition that 
the amplitude variation is very small, it does not matter much whether the 
bounce frequency is computed from the initial amplitude or from the time- 
asymptotrc amplitude. In our simulation, however, the amplitude variation 
is not negligible; if the bounce frequency were computed from the time-, 
asymptotic amplitudes, the points in Fig. 9 would be moved towards the 
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theoretical line of Morales and O'Neil. A similar amplitude dependence of 
the nonlinear frequency shift has recently been observed using the Vlasov 
approach. ^ 

Lee and Pocobelli have reported further simulation results of non- 
linear frequency shift.^ They used a one-dimensional particle simulation 

model (charge sheet model) with a spatial grid, and obtained agreement be- 

pS 

tween measured frequency shifts and theoretical calculations. In this 
simulation, a standing wave was driven by applying external electric field 
for a short time. This is in contrast to our case of an initial value prob- 
lem with a propagating wave. It is noteworthy that they used as many as 
40,000 particles. We needed only 4096 particles to obtain quantitative 
results of similar quality using the hybrid model. 
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6 . DISCUSSION 


In this part of our paper, we have studied essentially monochromatic 
electron plasma wave propagation by means of a variant of the economical, 
low-noise, hybrid simulation approach of Denavit.^" In Sections 3 and 4, 
computational studies of the model were made for equilibrium conditions, 
and small-signal propagation. The linear wave dispersion characteristics 
predicted by theory for long wavelength collective behavior were demon- 
strated to be followed very precisely. This verification of the validity 
and effectiveness of the simulation model is very important as a starting 
point for the subsequent study of nonlinear phenomena. It also serves to 

O 

establish the validity of the widely-used CIC model, and the Langdon 
theory describing the finite-size particle model. J Quantitative results 
in the very low energy range discussed here have never been obtained 
previously with such a high degree of accuracy with simple particle 
models . 


In the study of amplitude oscillation and Landau damping in Section 5, 
we have attempted investigation in areas where analytical approaches are 
not easily tractable, i.e., in cases where the condition, co/o^ » ( V p/ V t ) £ % 
is not satisfied. The results of our simulations show good qualitative 

agreement with the theories of Bailey and Denavit, 1 ^ and Sugihara and 

17 o 

Kamimura, who have made the assumption, oi/o^ » (v / v t ) t ' . However, 

there are significant differences between our simulation results and the 


theoretical results of these authors; first, phase-mixing of the amplitude 
oscillation is slower than predicted, and second, there exists no critical 
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value of 7^/0^ within our parameter range such as was found by Sugihara 
and Kamimura. It is hoped that these results will be helpful in better 
understanding the phenomenon , and in developing an analytical theory for 

s (v p/V 2 * 

In the study of nonlinear frequency shift made in Section 5, we have 
measured the frequency shift for large amplitude waves, and compared the 
results with theoretical predictions. It has been demonstrated that the 
simulation results agree well with the theoretical predictions of Morales 
and O’Neil. 27 

The deformation of the velocity distribution, particularly particle 
trapping, caused by a large-amplitude wave gives rise to a variety of wave 
amplification and coupling phenomena. Sidebands develop in addition to the 
monochromatic wave applied, and additional satellite frequencies appear. 
These can be studied as growth from noise, or by injection of suitable 
test waves. A series of such studies is presented in Part II of the 

35 

paper. 


3b 



REFERENCES 


1. J. Denavit, J. Comp. Phys. % 75 (1972). 

2. C. K. Birdsall and D. Fuss, J. Comp. Phys. 1+94 (19^9)* 

5- J. A. Byers, Proc. Fourth Conference on Numerical Simulation , 

November 1970 . Naval Research Laboratory, Washington, D. C., p. k$6 . 

4. J. M. Dawson, Phys. Rev. 118, 38 I (i 960 ). 

5 . H. R. Lewis, in Methods in Computational Physics , edited by 

B. Alder, S. Fernbach, and M. Rotenberg (Academic Press, New York, 

N. Y., 1970), Vol. 9 , p. 307 . 

6 . A. B. Langdon, J. Comp. Phys. 6 , 247 (1970). 

7- R. L. Morse and C. W. Nielson, Phys. Fluids 12 , 2418 (I 969 ). 

8 . T. P. Armstrong and C. W. Nielson, Phys. Fluids 1^, 1880 (I 97 O). 

9- J- Denavit and W. L. Kruer, Phys. Fluids _l4, 1782 (1971). 

10. Y. Matsuda, Stanford University Institute for Plasma Research Report 
No. 567 (March 1974). 

11. J. Canosa and J. Gazdag, Proc. Sixth Conference on Numerical Simulation 
of Plasmas. July 1973 . Lawrence Livermore Laboratory Report No. CONF- 
730804, Berkeley, California, p. 104. 

12. H. R. Lewis, Phys. Fluids Pj, 103 (1972). 

13 . T. M. O'Neil, Phys. Fluids 8 , 2255 (1965). 

14. L. M. Al'tshul and V. I. Karpman, Sov. Phys. JETP 22, 36 I ( 1966 ). 

15 . V. L. Bailey and J. Denavit, Phys. Fluids 1^, 451 (1970). 

16. S. P. Gary, Phys. Fluids 10, 570 ( 1967 ). 

17. R. Sugihara and T. Kamimura, J. Phys. Soc. Japan ^ 2 , 206 (I 972 ). 

18. I. H. Oei and D. G. Swanson, Phys. Fluids 1^, 2218 (I 972 ). 


35 



19- J. H. Malmberg and C. B. Wharton, Phys. Rev. Letters 12, 775 (I 967 ). 

20. A. Lee and G. Schmidt, Phys. Fluids 25 46 (1970). 

21. R. N. Franklin, S. M. Hamberger, and G. J. Smith, Phys. Rev. Letters 

22 , 91^ (1972). 

22. W. L. Kruer, J. M. Dawson, and R. N. Sudan, Phys. Rev. Letters 27 , 

838 (1969). 

23. G. Knorr, Z. Naturforsch. l 8 a , 1304 ( 1963 ). 

2k. T. P. Armstrong, Phys. Fluids 10, I 269 (I 967 ). 

25- J. M. Dawson and R. Shanny, Phys. Fluids ^L, I 5 O 6 (I 968 ). 

26. W. M. Manheimer and R. W. Flynn, Phys, Fluids _14, 2393 (1971). 

27. G. J. Morales and T. M. O'Neil, Phys. Rev. Letters 28, 417 (1972). 

28. R. L. Dewar, Phys. Fluids 1^_ } 712 (1972). 

29. A. Lee and G. Pocobelli, Phys. Fluids l^ } 2351 (1972). 

30. I. B. Bernstein, J. M. Green, and M. D. Kruskal, Phys. Rev. 108, 

546 (1957). 

31. W. E. Drummond and D. Pines, Nucl. Fusion Suppl . Pt. 3 , 1049 (I 962 ). 

32. A. A. Vedenov, E. P. Velikhov, and R. Z. Sagdeev, Nucl. Fusion Suppl. 

Pt. 2, 465 (1962). 

33* 0 • Canosa, private communication. 

34. A. Lee and G. Pocobelli, Phys. Fluids 16, 1964 (1973). 

35* Y. Matsuda and F. W. Crawford, Phys. Fluids [- Part II of the paper]. 


36 



